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We analyze the gapped phase of the Kitaev honeycomb model perturbatively in the isolated-dimer 
limit. Our analysis is based on the continuous unitary transformations method which allows one to 
compute the spectrum as well as matrix elements of operators between eigenstates, at high order. 
The starting point of our study consists in an exact mapping of the original honeycomb spin system 
onto a square-lattice model involving an effective spin and a hardcore boson. We then derive the 
low-energy effective Hamiltonian up to order 10 which is found to describe an interacting-anyon 
system, contrary to the order 4 result which predicts a free theory. These results give the ground- 
state energy in any vortex sector and thus also the vortex gap, which is relevant for experiments. 
Furthermore, we show that the elementary excitations are emerging free fermions composed of a 
hardcore boson with an attached spin- and phase- operator string. We also focus on observables 
and compute, in particular, the spin-spin correlation functions. We show that they admit a multi- 
plaquette expansion that we derive up to order 6. Finally, we study the creation and manipulation 
of anyons with local operators, show that they also create fermions, and discuss the relevance of our 
findings for experiments in optical lattices. 
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I. INTRODUCTION 

Elementary particles can be classified in two categories 
according to the value of their spin. Half-integer spin par- 
ticles obey Fermi-Dirac statistics and are called fermions 
whereas integer-spin particles obey Bose-Einstein statis- 
tics and are known as bosons. However, some quantum 
objects may obey other (fractional) statistics describing 
nontrivial braiding as initially suggested by Leinaas and 
Myrheim more than thirty years ago 1 , and by Wilczek 
in the eighties - ' ! . Despite numerous theoretical works, 
these so-called anyons are still waiting for a direct obser- 
vation although recent experimental proposals are very 
promising (see Ref. 4). 

In the last years, anyons have drawn much atten- 
tion because of their interest for topological quantum 
computation . In this perspective, several models have 
been proposed among which the celebrated toric code 1 ' 
which is a spin- 1/2 system whose elementary excitations 
behave as semions. However, the experimental realiza- 
tion of this system is rather tricky since it involves four- 
spin interactions. Here, we shall focus on another system 
originally proposed by Kitaev' which only involves two- 
spin interactions. This model is very rich since it contains 
Abelian and non-Abelian anyonic as well as fermionic ex- 
citations. Thus, it has been the subject of many recent 
studies concerning the spectrum, 8 ' 9 ' 10 ' 11 ' 12 ' 13 ' 14 ' 15 ' 16 the 
correlation functions and the entanglement 16 ' 17 ' 18 ' 19 ' 20 , 
or the quench dynamics 21,22 . Let us also mention sev- 
eral extensions 23,24 ' 25 among which the analysis of time- 
reversal symmetry breaking terms ' which may give 
rise to a chiral spin liquid. 

Furthermore, this model is susceptible to be realized in 
various experimental systems such as polar molecules, ul- 



tracold atoms 28 ' 29 ' 30 ' 31 or Josephson junctions 32 . It thus 
constitutes an appealing candidate for the observation 
of anyons. Nevertheless, the presence of fermions in the 
spectrum may spoil the detection process ; a point com- 
pletely missed in a recent proposal (see Ref. 33 for expla- 
nation and Ref. 16 for details). 

The goal of the present paper is to investigate the 
gapped phase of the Kitaev honeycomb model' . Indeed, 
in his remarkable seminal paper, Kitaev mainly focuses 
on the special subspace of the Hilbert space to which 
the ground state belongs to and the low-energy spec- 
trum of other subspaces has only been discussed lately 13 . 
Our aim is to bridge this gap by providing a high-order 
perturbative analysis, in the isolated-dimer limit, of the 
spectrum as well as some interesting results about the 
creation and the manipulation of anyons which is of rel- 
evance for experiments 30 ' 31 . Part of our results have al- 
ready been given in two short papers 13 ' 16 and the present 
paper may be considered as an extended and detailed ver- 
sion of these works. However many other results are pre- 
sented here among which the interplay between fermions 
and anyons under string operations discussed in Sec. IX. 

This paper is organized as follows. In the next section, 
we introduce the model as well as its main properties. 
In particular, we discuss the importance of the bound- 
ary conditions and insist on the role played by conserved 
quantities 1 " and the constraints resulting from them. In 
Sec. Ill, we show how to map the Kitaev model involv- 
ing spins on the honeycomb lattice onto an effective spin 
and hardcore boson on a square lattice. This mapping is 
the starting point of the perturbation theory presented 
in this work. In Sec. IV, we explain how to diagonal- 
ize the Hamiltonian order by order using the perturba- 
tive continuous unitary transformation (PCUT) method 
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. The study of the low-energy (zero-quasiparticle) sector 
is the subject of Sec. V, a large part of which is devoted 
to a pictural (and hopefully pedagogical) analysis and 
construction of the eigenstates of the toric code model 
which naturally emerges from this problem. There, we 
give the perturbative expansion form of the ground state 
energy for any vortex configuration. The effective low- 
energy theory is found to be described by interacting 
anyons contrary to the lowest-order result which predicts 
free anyons'. Sec. VI focuses on the study of the one- 
quasiparticle subspace, where the physics is shown to be 
that of a particle hopping in a magnetic field with zero 
or half a flux quantum per elementary plaquette. The 
demonstration of the fermionic nature (known from ex- 
act solutions) of the quasiparticles is briefly sketched in 
Sec. VII. In Sec. VIII, we provide some checks of our 
results by analyzing simple vortex configurations which 
allow for an exact solution. The spin-spin correlation 
functions and the manipulation of anyons are tackled 
in Sec. IX, which is devoted to the renormalization of 
observables. Finally, wc discuss several issues and give 
some perspectives. Technical details as well as all rele- 
vant coefficients involved in the perturbative expansions 
are gathered in appendices. 

In what follows, we tried to be as pedagogical as pos- 
sible and always favored simple demonstrations on con- 
crete examples rather than lengthy proofs for general sit- 
uations. We hope that it will help the reader to under- 
stand the richness of this model. 



II. THE MODEL 

A. Hamiltonian and boundary conditions 

The model considered in this work is a spin- 1/2 system 
proposed by Kitaev' in which spins are located at the 
vertices of a honeycomb lattice. Since the honeycomb 
lattice is topologically equivalent to the brick-wall lattice, 
we shall always represent it as shown in Fig. la. In this 
lattice, one distinguishes three types of links (x, y, and 
z) to which one associates three different couplings and 
interactions. The Hamiltonian of the system is 

a=x,y,z a — links 

where af are the usual Pauli matrices at site i. In the 
following we assume, without loss of generality , that 
J a > for all a and J z > J x , J y . 

We will cither work with an infinite system and open 
boundary conditions (a plane), or with a finite (or infi- 
nite) system and periodic boundary conditions (a torus). 
In the latter case and for reasons that will become clearer 
in the following (in particular, see Sec. VB), we shall 
restrict ourselves to the periodic boundary conditions 
(PBC) depicted in Fig. 2. The number of sites N s is 
N s = 2(2p) 2 = 8p 2 , with p e N (p = 1 in Fig. 2a). Let us 
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FIG. 1: (color online). Mapping of the honeycomb (brick- 
wall) lattice onto a square lattice with unit basis vectors rii 
and ri2. Each z— dirtier with four spin configurations is re- 
placed by a single site with four degrees of freedom : the 
occupation number of hardcore boson (0 or 1), and the effec- 
tive spin (ff, or ij.) which is chosen as the spin of the black site 
of the considered z — dimer. The numbering of the sites of a 
plaquette p is shown in both cases. 




FIG. 2: (color online). The periodic boundary conditions used 
in this work, (a) on the original brick-wall lattice and (b) on 
the effective square lattice of -2-dimers (see Fig. 1). In the 
figure p — 1. In both cases, the finite-size system is put on a 
torus obtained by identifying the opposite sides of the dashed 
(magenta) square. For clarity, the site at the point chosen as 
the origin has been depicted bigger (and in magenta). Half 
the square plaquettes in (b) have been colored in cyan (gray) 
to show that the periodic boundary conditions allow us to 
bi-color the lattice. 

anticipate what follows and mention that these boundary 
conditions are such that the lattice of z-dimers (Fig. lb) 
can be bi-colored as shown in Fig. 2b). 

B. Conserved quantities 

A remarkable property of Hamiltonian (1), is that its 
elementary operators Kij = erf a" commute with plaque- 
tte operators W p so that [H, W p ] — 0. For the plaquette 
p shown in Fig. la, such an operator is defined as 

W p = K 12 K 23 K 3i K i5 K 56 K 61 = <jla\al<rla\al. (2) 

Let us mention that in the expression of W p in terms 
of the K'b, one could have started at any site instead 
of site 1 and/or one could have taken the product of 
K's anti-clockwise instead of clockwise. Furthermore, 
the expression in terms of cr's could also be written W p = 
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FIG. 3: (color online). Illustration of the conserved plaquette 
quantity W p . The thick yellow (lightest gray) line delimitates 
the plaquette p. The thick red (gray), green (light gray) or 
blue (dark gray) segments represent the Pauli matrices a ° ut ( 1 ^ . 



out(i) 



where i runs over the set of six spins around 



the plaquette p, and where the notation out(i) means 
the "outgoing" direction at site i, with respect to the 
plaquette's contour. An illustration of the W p operator 
is given in Fig. 3. 

Since Wp = I, the eigenvalues of the plaquette oper- 
ators are w p = ±1. Note that [W p ,Wp'] = 0, as can 
be shown from the usual Pauli matrices algebra. As a 
consequence, H and the W p 's can be diagonalized simul- 
taneously. Following Kitaev, we will call a vortex sector 
a subspace of the Hilbert space with a given map of the 
Wp's. By definition a vortex is a plaquette for which 
w p = —1, so that for example, the vortex- free sector is 
defined by w p = +1 for all p's. 

In fact, all loop operators made of "outgoing spins" 
(see Figs. 4 and 5) are conserved and all commute with 
each other, which can be verified in the same way as for 
the Wp's. However, not all of them can be set indepen- 
dently to ±1. Some relations among them arise from 
the following fact (which can be checked by studying all 
possible cases) : the product of W p and a nearby loop 
operator C gives a new loop operator £ = W P C, as illus- 
trated on a particular example in Fig. 4. 

As an illustration of other relations involving loop op- 
erators around the torus, with the loops of Fig. 5 one 
has 
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FIG. 4: (color online). Illustration of the relation £' = W V C, 
with C = Yliec cr° ut ^ , W p already shown in Fig. 3 and 
C — Yliec' a° . The thick yellow (lightest gray) line in 
(a) represents the contour C and the one in (b) represents C' . 
As in Fig. 3, the thick red (gray), green (light gray) or blue 
(dark gray) segments represent the Pauli matrices CT ° ut W 
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FIG. 5: (color online). Examples of conserved loop operators 
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C = Yliac <r° ut ^ in a finite-size system with PBC. 



n=l 

8 

C' h =£ b l[W bn , 

n=l 

where we have denoted, for example, £ a = Yiiec CT 



(3) 

(4) 
(5) 

out(i) 



The minus sign in the last equation above comes from the 
crossing of £b an d C c . In the three expressions above, 
the product of plaquette operators could also have been 



taken over the complementary set of plaquettes. Indeed, 
on the torus the relations among loop operators yield the 
following constraint 



II w p 



I. 



(6) 



all p'i 



showing in particular that the number of vortices has to 
be even in a system with PBC. 

From examples shown in Fig. 5, one can deduce that 
all Wp's [except one, because of Eq. (6)], £b and C c can 
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be set independently to ±1, which then imposes all other 
conserved quantities. 

C. Some results from the exact solution 

The above discussed local conserved quantities are 
not sufficient to fully diagonalize the Hamiltonian. In- 
deed, if N s is the number of sites, then there is a total 
of TV — N s /2 plaquettes, but only N — 1 independent 
ones. With the two cycles around the torus, this gives 
N + 1 independent conserved quantities, which is obvi- 
ously smaller than N s . 

However, H has a crucial property : it can be trans- 
formed into a free Majorana fermions Hamiltonian and 
is thus exactly solvable. Let us also mention that an- 
other solution based on the Jordan- Wigner transfor- 
mation maps the spin Hamiltonian H onto a spinless 
fermions system with p-wave pairing 9 ' 11 ' 12 . 

As shown by Kitaev' , the ground state of H lies in the 
vortex- free sector and the phase diagram contains, in this 
sector, two phases : a gapped phase for J z > J x + J y and 
a gaplcss phase for J z < J x + J y . In the gapped phase 
the low-energy excitations are Abelian anyons (semions) 
whereas in the gapless phase, the low-energy excitations 
are fermionic. The gaplcss phase acquires a gap in the 
presence of a magnetic field and then contains gapped 
non- Abelian anyon excitations. The phase diagram has 
also been investigated in other vortex configurations such 
as the vortex-full sector and similar phases have been 
obtained. More precisely, one has a gapped phase for 
J z > J x + Jy and a gapless phase in the opposite case 14 . 

Our goal here is to determine the low-energy spectrum 
for any vortex configuration. Of course, one may use the 
fermionic Hamiltonian mentioned above but, it can only 
be exactly diagonalizcd for translation-invariant configu- 
ration. Here, we follow an alternative route by focusing 
on the isolated-dimer limit J z 3> J x , Jy 

III. MAPPING ONTO AN EFFECTIVE SPIN 
BOSON PROBLEM 

A. Mapping of the Hamiltonian 

The very first step of our analysis consists in mapping 
the four possible states of the two spins of a z-dimer onto 
those of an effective spin and a hardcore boson. More 
precisely, denoting | f) (| 1)) the eigenstate of a z with 
eigenvalue +1 (—1), an isolated z-dimer can be in one of 
the two low-energy states {| tT), 1 14)} with energy — J Zl 
or in one of the two high-energy states {| Tl)i I IT)} with 
energy + J z . Keeping in mind that our aim is to perform a 
perturbation theory in the limit J z 3> J x , J yi it is natural 
to interpret the change from a ferromagnetic to an anti- 
ferromagnetic configuration as the creation of a particle 
with an energy cost 2 J z . By construction, such a particle 
is a hardcore boson. The remaining degree of freedom 



can be described by a spin 1/2, indicating which of the 
two configurations is realized. There are many possible 
parametrizations but here we choose the following 

I TT> = |lro> ; |ii> = |4Lo>, |ti) = |fri),UT) = |4Li). (7) 

The left (right) spin is the one of the black (white) site of 
the dimer (| jj.) = | t»|o), etc). Double arrows represent 
the state of the effective spin which is here the same as 
the state of the left (black) spin. 

Within such a mapping, effective spins and hardcore 
bosons live on the effective square lattice of z-dimers (see 
Fig. 1). This lattice is shown again in Fig. 2b, together 
with the PBC, which are such that it can be bi-colored. 
In what follows, the sites of the effective lattice will be 
denoted with bold letters, such as i. 

Let us now write the Hamiltonian (1) in this language. 
Therefore, we first translate the action of the spin oper- 
ators in the effective-spin boson (ESB) formalism. It is 
easy to check that one has 

^l^rHbt + b,) , al^irfibt-b,), (8) 
°i,.=T? , ^,0=^(1-2^). 

The operators r" (a = x, y, z) are the Pauli matrices 
acting on the effective spin at site i, while b i and b\ 
arc hardcore bosonic annihilation and creation opera- 
tors, satisfying the usual on-site anticommutation rela- 
tion {foj,^} = I (and operators on different sites com- 
mute). Setting once for all J z = 1/2 so that creating a 
boson costs an energy 1 in the isolated-dimer limit, the 
Hamiltonian (1) reads 

H = -— + g + T +T +2 + T_ 2 , (9) 

where N is the number of z-dimers (or, equivalcntly, of 
square plaquettes), and 

Q = E 5 1^ ( 10 ) 

i 

T = - £ (j X + Jy + h.C.) , (11) 

i 

T + 2 = -E(^^ +Tll + J ^'' +ri2 ) = (T - 2)t - (12) 

i 

These operators arc built from local hopping and pair 
creation operators 

t i+ni = ht b-T? t i+n ' 2 = -\b ] - b-T? t? 

v? ni =4 +ni btT? +ni , v^=ib\ +n b\rV +n rl. 

(13) 

We emphasize that the mapping (8) explicitly breaks the 
symmetry between white and black sites of the original 
brick-wall lattice. This is responsible for the apparent 
breaking of symmetry between the xjri\ and yjni direc- 
tions in Eq. (13). However, for all the physically observ- 
able results, this symmetry remains intact (see for exam- 
ple the series expansion of eigenenergies in Appendix C). 
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Note however that the ri\ + rii and rt\ — ri2 directions 
are not equivalent, as can be seen from the underlying 
brick- wall lattice. 



B. Conserved quantities 

Let us now rephrase the conserved operators discussed 
in Sec. II B in the effective language. Using the nota- 
tions depicted in Fig. lb, as well as the mapping (8), the 
plaquette operators transform into 



W n = f-l) i, L i 'L+ 6 D b D T V T Z T V T Z 
rv P \ ) L U R D * 



(14) 



Note that (-l) b l b i = 1 - 2&tfe.. In the same vein, for the 
cycles around the torus shown in Figs. 5b-5c, which are 
reproduced for the effective lattice in Figs. 6b-6c, one has 



A> = n, 



-(-1) 



6Tb. 



n 



there is an even number of sites on the contour with 
the PBC chosen here), as well as £ c = Yliec T t ■ ^ nc 
expression for Cd (see Fig. 5d), namely Cd — IlieCd^' 
is a bit more complicated, but it should be clear from 
Fig. 6d what the uVs are. Finally, for the contour shown 
in Fig. 6a (which is in correspondence with Fig. 5a), one 



has Y 



pcc t 



n 



i6C a 



= C a with the Wj's indicated 



in the figure, and with p C C a meaning the plaqucttes 
p enclosed in the contour C a . With these notations, one 
can easily check that Eq. (5) still holds. 

The elementary hopping and pair creation operators, 
namely t\ and v\ with i and j nearest neighbors, have 
a very remarkable property : they all commute with the 



Wp's, as well as with any other loop operator 



[1Z,W p ] = [v>,W p ] 



o. 



(15) 



The original spin problem on the honeycomb lattice is 
thus mapped onto a quadratic hardcore boson problem 
on an effective square lattice, with conserved plaquette 
and loop operators. Let us underline that this mapping 
is exact and just provides an alternative description of 
the spin problem. The resulting Hamiltonian (9) remains 
difficult to diagonalizc (except, of course, if one remem- 
bers that the model can be fermionized), since (i) bosons 
are hard core which prevents the use of a Bogoliubov 
transformation, (ii) bosonic and spin degrees of freedom 
are correlated. The conserved plaquette operators will of 
course be useful in simplifying and solving the problem 
as recently underlined in Ref. 15. 



IV. PERTURBATION THEORY IN THE 
GAPPED PHASE 

A. Effective Hamiltonian from PCUTs 

The starting point of the present perturbation theory 
is the isolated-dimcr limit, namely J x = J y = 0. In this 







FIG. 6: (color online). On the four figures, the thick yellow 
(lightest gray) line represents the contours C. The operators 
uii which are such that the loop operators read C = Hiec^ 
are indicated in the figures. In the present case, they can take 
the following values x = r x , x — ( — l) b b r x (and the same for 

y and z), and —x — — ( — l) 6 b r x . These figures are the same 
as the ones of Fig. 5, but on the effective lattice. 



limit, the spectrum is made of equidistant and degenerate 
levels separated by an energy gap A = 2J Z = 1. To com- 
pute the perturbative spectrum, there are of course sev- 
eral methods among which the Green's function formal- 
ism initially used by Kitaev' . However, if this approach is 
efficient to obtain the first nontrivial (nonconstant) cor- 
rection, it becomes tricky to implement at higher orders. 

Here, following Ref. 13, we use an alternative ap- 
proach based on continuous unitary transformations 
(CUTs) conjointly proposed by Wegner 34 and Glazek and 
Wilson 35 ' . We refer the interested reader to Ref. 37 
for a recent pedagogical introduction. Its perturbative 
version denoted PCUTs is especially well-suited to the 
problem at hand. This technique is detailed in sev- 
eral works 38 ' 39 . Let us simply mention that the CUTs 
method requires the choice of a generator that drives the 
flow of the operators. All the results given here have 
been obtained with the so-called quasiparticle number- 
conserving generator first proposed by Mielke 40 for finite 
matrices and generalized to many-body systems by Knct- 
ter and Uhrig 39 . 

The latter have computed the perturbative expansion 
for any Hamiltonian of the form 



T- 



-2- 



(16) 



H = Q + T-2 + T_i + To + T +1 
provided two hypothesis are satisfied: 

• the unperturbed Hamiltonian Q has an equidistant 



6 



spectrum bounded from below ; 

• the perturbing Hamiltonian ^2n=-2 -^n is such that 
[Q,T n ] = nT n . 

Clearly, the Hamiltonian (9) meets these two criteria (up 
to a constant term) noting that in the present case, one 
has T±i = 0. Here, we have included the "small" pa- 
rameters, namely J x and J y , in the definition of the T n 
operators, which is not the convention usually adopted 
in the CUTs community. 

The CUTs method together with the quasiparticle 
number-conserving generator unitarily transforms the 
Hamiltonian (16) into an effective Hamiltonian H e s = 
U* HU commuting with Q, U being a unitary operator. 
We give the first terms of the expansion up to order 4 
in Appendix A. As can be seen in Table I, the number 
of terms appearing in the perturbative expansion quickly 
increases with the order. For instance, at order 2, the 
effective Hamiltonian reads in our case 

N 11 

H eS = -— + Q + T - -T_ 2 T +2 + -T +2 T_ 2 , (17) 

whereas at order 10, there are more than 10 4 operators 
to consider. 

Writing H e s this way is only the very first part of the 
job since one next has to (i) determine its action in each 
subspacc of a given quasiparticle (QP) number q ; (ii) 
diagonalize H c ff in each of these subspaces. This is the 
object of the next sections : we first study the lowest- 
energy states (q = QP) which is the main contribution 
of our work ; then we turn to the q = 1 QP states and 
recover the high-energy gap from the QP dispersion ; we 
end by q ^ 2 QP states, whose properties determine the 
statistics of the QPs, and we will see that the QPs behave 
as fermions which arc, furthermore, non interacting. This 
fact is at the origin of a tremendous simplification of the 
effective Hamiltonian. Indeed, we found that H e ff can be 
written, at all orders and in the thermodynamical limit, 
as 

H cS = E +nQ- J2 C vu-,PnW Pl . . .W Pn (18) 

{Pl,—,Pn} 

E i h J h :ir 

{ii>—.i„} 

We shall discuss each term in detail in the following 
sections, but let us mention that Eg, /i, the C"s and 
the D's are coefficients whose series expansion are com- 
puted. The S operators are strings of spin operators rj* 

and of phase factors (— l) b j b i on the cluster . . . , j n }. 
This very special form of multi-particle terms [remember 

(— l) b 3 b j = 1 — 2£>t&.], leading to phase factors and spin- 
strings only, is responsible for the emergence of fermions 
in the model. 

For a finite-size system with PBC, new terms appear 
in the effective Hamiltonian. They involve loop operators 
around the torus, and appear at a minimal order being 



the linear size 2p of the lattice. Such loop operators are 
associated to contours as the ones shown in Figs. 5 and 
6, namely £b, £c and C' h for the contours Cb, C c and 
C b . The presence of such loop operators in the effective 
Hamiltonian shows that the eigenstates of the Hamilto- 
nian are also eigenstates of these loop operators. Their 
effect is to lift the degeneracies between states (which 
for each energy is at least four in the thermodynamical 
limit, since some of the excitations are Abelian semions 
and the genus of a torus is 1, see Ref. 5). We shall not 
dive into the details of such finite-size corrections, since 
our approach allows us to directly tackle with the most 
interesting thermodynamical limit. However, let us make 
a remark about a numerical check of this statement for 
small system sizes. For a torus whose linear size is strictly 
smaller than 4, the loop operator terms around the torus 
dominate the expansion over the W p 's, and for a size of 4 
both types of terms start contributing at the same order. 
One should thus not be surprised to find a ground-state 
for p = 1 which is not in the vortex- free sector . 



B. Counting of states 

Before we turn to a detailed analysis of each QP sub- 
space, let us show that we do not miss any state us- 
ing simple counting arguments. We have already seen 
in Sec. II, that one has A + 2 conserved Z2 quantities 
(two loop operators, and N plaquette operators), with 
the constraint rj W p = I. There is in fact one more re- 
lation between the W p 's, involving the number of bosons, 
which reads 

n w P ={-i)^ b ^={-\f= n w P , 

pewhite pecyan (gray) 

(19) 

showing that the parity of the number of vortices living 
on white plaquettes (see Fig. 2) has to be the same as the 
parity of the number of bosons. The last equality simply 
comes from the previously mentioned constraint (6). The 
first equality can be checked using the expression (14) of 
the Wp's. Indeed, for a site i having a white plaquette on 
its left, and another one on its right, the product of the 
two associated W p 's will give if x (-l) 6 * b *7f = (-l) fc I b ;. 
In the same way, for a site i having a white plaquette 
above it, and another one under it, the product of the 
two associated W p 's will give r? x (— lj^^T? = (— \) b ^ b i . 
Let us note that (19) has a meaning in the two bases we 
are working in, the initial one and the unitarily trans- 
formed one. Indeed, in the initial basis, the Hamiltonian 
H commutes with the parity operator (— 1)^ ; in the 
rotated basis, H e g commutes with Q. 

We thus see that in a subspace with a given number of 
QP's, there are N independent conserved 1 2 quantities. 
Thus, N being the number of effective spins Ti, there is 
no remaining effective spin degree of freedom once the 
1i 2 quantities are chosen. As a conclusion, the g-QP sub- 



7 



space has dimension d q = 2 ( g J ' w ^ n the usual no- 
tation for binomial coefficients. This shows that we miss 
no state, since 

N N / \ 

E d * = 2W Eu)= 22Ar = 2JVB > ( 2 °) 

where N s is the total number of spins in the brick-wall 
lattice. 

This discussion furthermore sheds light on the fact that 
in order to compute eigenenergies, a perturbative expan- 
sion of the Kitaev model (as opposed to exact numerics) 
is really of interest only in the 0-QP subspace. Indeed, 
we have just seen that there are N independent Z 2 con- 
served quantities. It is thus clear that as soon as we 
will have written down the effective Hamiltonian in the 
0-QP subspace, the Hamiltonian will already be diago- 
nal, whatever the vortex configuration, although writing 
down the eigenstates of the Z 2 quantities in the basis of 
effective spin operators still has to be done. However, in 
the 1-QP subspace, one will have to diagonalize an N x N 
matrix (numerically in the case of a nonperiodic vortex 
configuration), which is identical to what one has to do 
when solving the problem exactly as Kitaev did. For 
q ^ 2, the perturbative expansion looks even more com- 
plicated than the exact solution, but this is an artifact, 
since we recover free fermions. 



V. EFFECTIVE HAMILTONIAN IN THE 0-QP 
SUBSPACE 

A. Effective Hamiltonian and eigenenergies 

In the 0-QP sector, and in the thcrmodynamical limit 
the effective Hamiltonian (18) simplifies and reads 

H cS \ q=0 = E - C Pi,-,PnW Pl ...W Pn , (21) 

{pi,...,p„} 

where {pi,p2, ■ ■ ■ ,Pn\ denotes a set of n plaquettes and 
the Wp's are the conserved plaquette operators intro- 
duced in Sec. II. Note than when restricted to the 0- 
QP sector they simplify to Wp| g=0 = t^t^t^t^ [see 
Eq. (14)]. 

As mentioned at the end of the previous section, ob- 
taining eigenenergies only requires a minimal amount of 
work, namely replacing each W p by numbers w p = ±1, 
and doing the same with loop operators, without for- 
getting about the constraints among these quantities. 
The perturbative expansion of the coefficients Eq and 
Cpi,...,p„ are given in Appendix C. Let us note that 
{pi,p 2 , • ■ • ,Pn} does not need to be a linked cluster of 
plaquettes (as seen for C PjP -\-2ni that is nonvanishing at 
order 10), and that translational invariancc of the Hamil- 
tonian implies that the C pi ^ ^ Pn coefficients only depend 
on n — 1 relative positions of the plaquettes. 
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FIG. 7: (color online). Two-anyon configurations (gray cen- 
tral plaquette and one of the numbered plaquettes) on a 
vortex-free background. A_Bi v (Ai?2v) is the energy cost (at 
order 10) for adding one vortex (two vortices) to the vortex- 
free state. A£iv reads A£i v = J 4 + 8J 6 + 75 J 8 + 784 J 10 . 
For simplicity, we have set here J x = J y = J but the results, 
in the general case, are easily obtained from the coefficients 
given in Appendix C. 



The lowest nontrivial order involving the Wp's (order 
4) has been derived by Kitaev' (C p = J^Jy/2) and led 
him to identify the effective low-energy theory with the 
toric code . One of the main results of our work is to 
show that, at order 6 and beyond, one obtains a multi- 
plaquette expansion in the effective low-energy Hamilto- 
nian. In other words vortices interact, though they re- 
main static as they have to since the Wp's are conserved. 
The interaction energies between vortices are not directly 
the C coefficients. One should write w p = 1 — 2n p where 
rip is the number of vortices at plaquette p (0 or 1), then 
look at coefficients in the expansion in terms of the n p s. 
The results of such an analysis for two-vortex interac- 
tion energies in the case J x = J y = J are illustrated in 
Fig. 7 which shows that the interaction (i) lowers the en- 
ergy and is therefore attractive, (ii) is anisotropic even 
for J x = J y = J which is clear from the structure of 
the underlying brick- wall lattice, (in) decreases with the 
distance d between vortices as expected in a gapped sys- 
tem. Note that for a finite-size system with PBC, the 
two-vortex configurations with a central vortex and an- 



FIG. 8: (color online). Illustration of the two different points 
of view one can have of a bi-colored lattice : sites are at 
vertices in (a) and on the bonds in (b). The plaquette m in 
(a) remains a plaquette in (b), while the plaquette operator 
e transforms into a star operator. 



other vortex at one of 1, 5, 6 and 7 sites (see Fig. 7) are 
forbidden since they violate the constraint (19). A one- 
vortex configuration is also forbidden since it violates the 
constraint (6). These configurations would be allowed in 
an infinite system, or in a finite system with open bound- 
ary conditions. 

A most remarkable point which emerges from the anal- 
ysis of i/eff |o is that its eigenstates are those of the Wp's. 
They are thus the same at any order (> 4), and are those 
of the toric code (> , although their eigenenergy changes 
with the perturbation order (we emphasize we are talking 
about eigenstates of H e g |o, not of the original Hamilto- 
nian H). We graphically sketch the construction of these 
eigenstates in the next subsection, which will also prove 
to be useful for the (q ^ 1) QP sectors, and show ex- 
plicitly that they obey anyonic, more precisely semionic, 
statistics. Our discussion of the toric code focuses on 
peculiarities related to our way of studying the problem, 
that is not restricted to the 0-QP subspace. For more de- 
tails about the toric code model, we refer the interested 
reader to Refs. 5,6,7. 



B. The Toric Code in a nutshell 

1. Mapping to the toric code 

As we have seen in the previous sections, the eigen- 
states of the effective Hamiltonian in the 0-QP subspace 
are the eigenstates of the Hp's and of the £'s. We re- 
call that in this subspace, the plaquette operators read 
W p \ q=0 = tI t* tI t* (sec Eq. (14) and Fig. lb). A sim- 
ilar simplification occurs for the £'s. As mentioned by 
Kitaev' (for the Hamiltonian at order 4), the effective 
Hamiltonian could be studied directly, but it is much eas- 
ier to visualize the eigenstates by performing some spin 
rotations, and bring the Hamiltonian to the one of the 
toric code (generalized by multi- vortex terms). Thanks 
to the special PBC we have chosen, the lattice sites can 
be bi-colored in black and white as illustrated in Fig. 8. 
Then, one performs a different rotation on the two kinds 



FIG. 9: (color online). Two of the states entering the equal- 
weight superposition in (24). Crosses on the sites indicate 
a spin flip, with respect to the reference state | "ff) where all 
spins point upwards (the A e operators are still pictured by 
thick crosses at the vertices). The loops, in dotted lines, are 
obtained by joining the flipped spins. 



of sites 




This way, a cyan (gray) (resp. white) plaquette such 
as m (resp. e) in Fig. 8a transforms into a plaque- 
tte (star) term B m = ( — l) fc L fc L+ & D fc D s * s * s * s * (resp. 
A e = (— l) & L h L+ b D fc D s * s x s £ s £) 5 as shown with thick 
(red) lines in Fig. 8b. We have kept track of the phases 
involving boson numbers because our construction will 
be needed for (q ^ l)-subspaces, but it is clear that they 
can be dropped in the 0-QP subspace. Let us mention 
that the distinction between plaquette and star terms is 
purely conventional. The letters m and e refer to the 
magnetic and electric vocabulary also used by Kitaev', 
although we emphasize there is absolutely no difference 
between an A and a B operator, which are both disguised 
W operators. Up to an additive constant term, the ef- 
fective Hamiltonian in the 0-QP subspace, and at order 
4 finally reads (with J off = J^J^/2) 

H cS \ q=0 = -Jeff I J2 A e+J2 Bm ) ■ ( 23 ) 
\ e m / 

We work with this lowest (nontrivial) order Hamilto- 
nian, because the eigenstates of iJ c ff |<j=o remain the same 
whatever the order in perturbation. One should simply 
remember that the eigenstates of H also have to be eigen- 
states of C operators. If the PBC are not of the type we 
use (see Fig. 2), the sites can usually not be bi-colored 
and the rotations (22) cannot be performed, which makes 
the construction much more complicated, and we shall 
refer the interested readers to Ref. 41. 



2. Construction of the ground state (s) 

As a warming up, let us construct a ground-state of 
J/efflcb e -j an eigenstate of all B m and A e operators 



FIG. 10: (color online). The two loop operators C x and C v 
and their contours in yellow (lightest gray) thick lines, corre- 
sponding to Cb and C c in Fig. 6. As in the latter figure, but 
for s spins instead of r spins, y means s v , etc. In the 0-QP 
sector, y = {—V) hb s y is the same as y. 

with eigenvalue 1 (there are actually four of these). An 
eigenstate of all i? m 's is for example the "reference" state 
| ft) where all spins s point in the +z-dircction, such that 
for all i one has sf | ft) = | ft) . This state is not an eigen- 
state of the A e operators yet, but a simple projection 
yields the desired state 

2* /4 - 1/2 n(^V>®l°>>> ( 24 ) 

whose normalization follows from the number N/2 of e's 
and the property ]J e A e = I [see Eq. (19)]. The state |0)h 
indicates that there is no quasiparticle, i. e. no hardccore 
boson. A graphical interpretation can be given of the 
state (24) : it is an equal-weight superposition of multi- 
loop configurations produced by the A e operators, as the 
ones shown in Fig. 9. 

One next has to get an eigenstate of two independent 
loop operators, which we choose to be Ch and C c (see 
Figs. 5 and 6) and which, from now on, will be denoted 
C x and C v , with eigenvalues l x and l y . The expressions 
of these operators in the s-spin language are given in 
Fig. 10. 

Note that in the 0-QP subspace, one could also have 
used other conserved loop operators which are products 
of s x or of s z on the contours defining C x and C v . Such 
operators resemble more the ones used by Kitaev (> , but 
they are conserved only in the 0-QP subspace (in con- 
trast to C x and £ y ), and so will not prove to be very 
useful in the following. As can be seen in Fig. 10, C x 
and C y perform spin flips, with respect to | ft) on their 
associated contours. The four ground-states of (23) are 
then obtained with another projection and proper nor- 
malization 

|{ Wp = l},P,F) = 2 w / 4+1 / 2 (25) 
fl + l x C x \ fI + l y C y \ tj fI + A e \ 1A . .. 

x [—^- j j n [— j i ft) ® io), 

These four states are equal-weight (in absolute value) su- 
perposition of all possible multi-loop configurations, pro- 
duced by the A e operators as in Fig. 9, as well as the C x 
and C v operators, as illustrated in Fig. 11 for C x . 



FIG. 11: (color online). Two of the states entering the equal- 
weight (in absolute value) superposition in (25), involving the 
C x loop operator. Graphical conventions are the same as in 
Figs. 9 and 10. 



Let us note that the preceding construction relies on 
the | ft) state and the fact that it is an eigenstate of the 
_B m 's, etc. However, one could also have started with a 
state | =>■) where all spins point in the x-direction, which 
is an eigenstate of the A e 's and then follow a similar- 
route. 



3. Construction of excited states 

We now have to see how to construct excited states, 
i. e. states containing vortices (e or m) but still no quasi- 
particle. 

Constructing a state with some 5 m 's being minus one 
("magnetic vortices") is easy, once one has noticed that 
s x anticommutes with two J3 m 's, and thus changes their 
values to their opposite. Since s x commutes with all A e 's, 
as well as with C x and C v when i does not belong to the 
corresponding contours, sf \{w p = 1}, l x , I v )q is an eigen- 
state of the effective Hamiltonian, with two vortices living 
on the plaquettes touching the bond to which i belongs. 
(Note that since s x also anticommutes with C x and C v 
when i belongs to the corresponding contours, one should 
use a string of s x going around the torus without cross- 
ing C x and C y instead of s x ). The corresponding state is 
again an equal-weight (in absolute value) superposition 
of states, but now with all possible open strings joining 
the created vortices, as well as all possible closed loops. 
This is illustrated in Fig. 12. 

Creating "electric" vortices is as easy, since one can 
replace EL (^) m Eq. (25) by EL ( I± %^), with a e = 
±1, respecting the constraint Ile a e = 1 [ see Eq. (19)]. 
Such a change can also be obtained via the action of sf 
operators. Indeed, each sf operator anticommutes with 
two AgS, and thus changes their values to their opposite. 
The fluctuation of the strings induced by sf operators 
is however hard to see with the construction we have 
given, which relies on the reference state | ft). To see 
this, one should construct states from the reference state 
| =>•) where all spins point in the x direction, and then 
use projectors involving -B m 's instead of yl e 's. This is 



10 




FIG. 12: (color online). Two of the states entering the equal- 
weight (in absolute value) superposition for the state having 
two "magnetic" vortices. The latter are represented with little 
gray squares marked with the letter m, and are linked with 
a string of spin- flips [blue (dark gray) thick line] . The other 
graphical conventions are the same as in Figs. 9. 



not useful for our purpose so we let the interested reader 
doing it on his own. 




FIG. 13: (color online). Illustration of operators involved in 
the braiding of a "magnetic" vortex around an electric vortex. 
X is the product JT^ sf for i belonging to the blue (dark gray) 
thick path linking the two m vortices (orthogonal to bonds). 
Z is the product sf for i belonging to the green (light gray) 
thick path linking the two e vortices (drawn on the bonds). 
X' is the product JT f sf for i belonging to the dashed loop 
around the leftmost e vortex, and is also equal to a product 
of the A e 's operators encircled by the loop (denoted as thick 
crosses) . 



4- The statistics of vortices 

For completeness, let us now show that "magnetic" and 
" electric" vortices behave as semions with respect to each 
other. This is done by first creating a pair of "magnetic" 
vortices, then a pair of "electric" vortices, and finally by 
moving one of the " magnetic" vortices around one of the 
"electric" vortices as shown in Fig. 13 (one could also do 
the contrary, but then one should work with the refer- 
ence state | =>) to see things more easily) . With the no- 
tations of this figure (see also its caption) , let us consider 
the state = ZX\{w p = l},l x ,l v )o with two e and 
m vortices. Then the repeated application of spin-flips 
along the loop X' (in any direction) moves the down- 
most m vortex around the leftmost c vortex. The result- 
ing state is X'\yb). But as Z and X' have one (and only 
one) common site, they anticommute, whereas X and X' 
commute, so that X'\ip) = -ZXX'\{w p = 1},F,^) . 
Now, X' which is a product of s| operators forming a 
closed loop is nothing but a product of A e 's operators 
(the ones enclosed in the loop). As \{w p = l},l x ,l v )o is 
an eigenstate of the A^s with eigenvalue one, we finally 
obtain that X'\tp) = —\ip) : braiding a magnetic vortex 
around an electric vortex yields a nontrivial phase of it 
(— 1 = e 1T ), which proves the semionic statistics. 

Let us mention that the magnetic vortices behave as 
bosons among themselves, and so do the electric vortices. 
This is easily seen by noticing that creating and moving 
m vortices for example, only requires s x operators, which 
all commute with one another. To end this discussion 
about the statistics of vortices, let us also remark that 
a compound object made of an electric and a magnetic 
vortex is a fermion (see Ref. 5). 



VI. EFFECTIVE HAMILTONIAN IN THE 1-QP 
SUBSPACE 

A. Form of the Hamiltonian 

The spectrum we obtained in the 0-QP subspace gives 
the lowest eigenenergies for each configuration of the 
Wp's. In this section, we explain how to compute the 
high-energy spectrum for states with one quasiparticle, 
for each W p 's configuration, and how to build the associ- 
ated eigenstatcs. This is achieved by diagonalizing H c g 
in the 1-QP subspace (whose dimension is d\ = N2 N , 
see the end of Sec. IV). In this subspace, the effective 
Hamiltonian (18) reads 

H eS \ q =i = E + pL- ^2 C pu-,Pn W Pi ■■■ W Pn 
{pi,...,p„} 

E D i J S :i j >'] >' :i ■ (26) 

tii i„} 

where the second sum is performed over all non self- 
retracing paths of length n starting at site j 1 and ending 
at site j n , with possibly j n — j 1 when working at order 
4 or higher. This is the reason why we give the expan- 
sion up to this order but we would like to emphasize 
that obtaining orders up to 10 for H c r \ q= \ is of the same 
complexity as for -ff c ff |<?=o- Self- retracing paths are renor- 
malizing the chemical potential [i. Note that a hopping 
process of one quasiparticle around a loop is nothing but 
the product of the Wp's enclosed in the loop, as can be 
easily checked. This explains why at order 4, one obtains 
some terms proportional to bJ>-W p , where the plaquette 
p shares site i (see Appendix D). 

From now on (q 1) the phase factors appearing in 
the Wp's [see Eq. (14)] must be taken into account. The 
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operators 5 have a structure similar to that of the W p 's, 
except that they are open string operators. They involve 

Tj as well as phase factors (— 1) j a as follows 

*j , w j V; ■■■V'- (27) 

where the are phase factors which reduce to the 

identity in the 1-QP subspace, and will be discussed later 
on (see Sec. VII). The two-site TP operators are built 
from the same rj* operators as the hoppings tP , namely 

V +ni = rt +ni = {T^, (28) 
V +n2 = -^ +n2 Tf = (V +n J. (29) 

Note that in the 1-QP subspace, one can also write the 
hopping term of the Hamiltonian as 1. 




FIG. 14: (color online). Illustration of states of the 1-QP 
basis, built from the ground-state. The action of is to 
create a particle at site O (large black filled circle) as well as 
to create a pair of e and m vortices (a). The action of the 
string operator S represented with an oriented thick (cyan) 
line, on the state (a), then yields the state with the same 
vortices, but a particle at site i (b). We do not bi-color the 
lattice in cyan (gray) and white any more, so that the figures 
are easier to read. 



B. Construction of a 1-QP basis 

As can be seen when looking at the form of the hop- 
ping operators, bosonic and spin degrees of freedom are 
coupled so that one has to tackle a polaron-like problem. 
However, we shall now show that since all hopping oper- 
ators tP commute with all W p s as well as with all loop 
operators C, the 1-QP problem is equivalent to that of 
one particle hopping in a static magnetic field. 

As a first step, we build a basis of the 1-QP sub- 
space. We denote by \{w p },l x ,l y }o a state of the 0- 
QP subspace, which is an eigenstate of the W p s and 
of C x and C v , and built as explained in Sec. VB. We 
choose as the origin O the site we have already denoted 
with a large (magenta) filled circle (see Figs. 2 and 6, 
as well as Fig. 14b). Let us then consider the state 
\{w' p },l x 1 iy-0) l = bl\{w p },l x ,iy) belonging to the 1- 
QP subspace, and with one QP at the origin. From for- 
mula (14), it is clear that adding a particle at the ori- 
gin changes the value of two plaquettes, as illustrated in 
Fig. 14a for the action of b Q on the ground state, which is 
the reason why we made a distinction between w p and w' . 
Note that all this is perfectly consistent with Eq. (19), 
as well as with the conclusion of Ref. 42. Indeed, in this 
paper, Levin and Wen showed that fermions are always 
created in pairs, and this is the case here since a bound 
object of an "electric" vortex and a "magnetic" vortex is 
a fcrmion (see Sec. VB) and our quasiparticles will turn 
out to be fermions (see Sec. VII). 

Other states \{w' p } , l x , l v ; i) i with a particle at site i 

arc obtained by applying an operator So ib\b a onto 
l x , l v ; 0}i in order to make the particle hop, with- 
out affecting the conserved Z 2 quantities. Note that 

So,....ib%\{w' p },l x ,iy-0) 1 = 5 ,...,4lK}> ^>o, 

(31) 



However, we still need a convention for the path to 
be taken (which will amount to choose a gauge for the 
magnetic field the particles arc hopping in) to obtain a 
well-defined basis. The path from O to i is taken to first 
be in the rt\ direction as much as needed, then in the n 2 
direction. For example, in Fig. 14b, i = 3ni + 2n.2 and 
the 5 operator is depicted as an oriented thick (cyan) 
line in this figure, with first 3 moves in direction ri\ and 
then 2 moves in direction n 2 . 



C. Hamiltonian in the 1-QP basis 

Let us now consider the effective Hamiltonian at 
order 1, for which the hopping part is nothing but 
T , and study its action on a state \{w p },l x ,l v ;i)i. 
From the way the states have been built, it is obvious 
that tf n ^\{w p },l x ,iy;i)x = \{w p }, l x , l v ; i + n 2 )i, and 
for the same reason and the fact that TP is unitary, 
tr n2 \{w P },l x ,l y ;i)i = \{w p },l x > l v ;i-n 2 )i. We then 
turn to the hopping term tl +ni and study its action on 
the state \{w p },l x ,l v ;i)i. In other words, we wish to 
compute the matrix element 

^• +rM = i({w p },l x : iy-i + ni \ti +ni !{«;„}, I*, I*; 

(32) 

All needed states are represented in Fig. 15 : 
\{w p },l x jy;i) 1 in (a), t^ ni \{w p }, l x , l v ; t)i in (b) and 
\{w p } , l x , l v ; i + ni)i in (c). Using the notations of 

Fig. 15 (5 a is the oriented product of spin operators T? 
on the contour shown in (a), starting at the origin, the 
same for 5b and 5 C , but for 5d the product starts and 
ends at the particle's position), it is easy to see that 
S d t i i +n i\{w p }J x ,iy;i) 1 = \{w p },l x ,l y ;i + n 1 ) 1 . Then, 
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FIG. 15: (color online). States (a), (b), (c) and contour (d) 
needed to compute the matrix element (32), see text. 

using the fact that = I, 

K +ni = i({w p },l x ,l v -,i + ni\S d \{w p },l x ,I a ;i + n 1 ) 1 , 

(33) 

Furthermore a calculation on Pauli matrices shows that 
the action of Sd on the state l x , l v ; i + is the 

same as the product of the plaquette operators encircled 
by the closed contour of which on the example of 
Fig. 15 reads W P1 W P2 W P3 . We finally obtain 

A^= [] ( 34 ) 
pc£i,i+ ni 

where the product has to be taken over all encircled pla- 
quettes £i,i+m as illustrated on a particular example in 
Fig. 15. The case of a hopping in the —n\ direction can 
of course be deduced from the above matrix elements by 
hermitian conjugation. 

For some hopping processes, the matrix element not 
only involves a product of w p s but also a loop operator 
around the torus. This is illustrated in Fig. 16 for a 
hopping in the n 2 direction, starting from site i = i x n\ + 
(2p — l)n 2 (i x — 2 and p = 2 in the figure). In this 
case, one has S<i — L\ (i. e., the loop operator in the 
y-direction, around the torus, going through the sites i = 
i x rii + pn 2 where p takes all possible values), so that 

Af n ' 2 = i (K>, i", i y -A tf n2 1 W, i x , i y ; i)i = % ■ 

(35) 

As explained in Sees. II and III the value of l\ is deter- 
mined from the one of l v = l\ and from the value of the 
plaquettes in between these two loop operators. 




FIG. 16: (color online). States (a), (b), (c) and contour (d) 
needed to compute the matrix element of a hopping in the n2 
direction, when site i is "one site away from the edge" of the 
dashed (magenta) square, and thus involving a loop operator 
(d). See text for details. 



All the above examples lead to the following conclu- 
sion. The matrix elements of the effective Hamiltonian 
at lowest order, in the TQP subspace, are the same as 
what one would obtain for a particle with hopping am- 
plitudes — J x and —J y in the n\ and n 2 direction of the 
square lattice, in a magnetic field whose (reduced) fluxes 
in plaquettes or cycles around the torus are <E» / <3>o = 
or Q/$o = 1/2 (where $o is the flux quantum). This 
comes from the fact that, for example, hopping around a 
plaquette p gives a phase factor 1 for the two hoppings in 
the ±ri2 directions, and an overall w p for the hoppings 
in the ±ni directions. The overall contribution is then 
w p , which takes value w p = exp[2i7r$/$o]- This analy- 
sis can be extended to the case of hoppings of the kind 
represented in Fig. 16 where the PBC play a role. 

When tackling higher-order corrections, hoppings be- 
come longer-ranged as seen in Eq. (26), but the above 
considerations still apply because of Eq. (30). It is then 
easy to compute the 1-QP spectrum for a given map of 
the Z2 conserved quantities. As already explained, when 
the map docs not possess translational invariancc, one 
can only compute the spectrum numerically. When the 
w p s are translationally invariant, an analytic solution is 
available, and for example, in the vortex-free subspace, 
the dispersion relation obtained at order 2 (see Appendix 
D) is 

E hcc {k Xl k y ) = l-2[J x cos{k x ) + J y cos{k y )] (36) 
+2[J X s'm(k x ) + J y sin(fcy)] 2 , 
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FIG. 17: (color online). Illustration of the phase factors ap- 
pearing at site j 2 in the string operators S...j 1 j 2 j 3 ,.... The 

phase factors ( — l) b i2 bj 2 are denoted as large (cyan) dots 
and are only involved in the three top processes which are, 



from left to right, S. 
and S. . 7 ,7 —no . 7 +m — n 



■,3,3 

I,--- 



-™1 ,J+2tii , 



S. 



■J ,3+™2,3 +n, 1 +n 2 ,. 



where the wave vector (k x ,k y ) belongs to [— tt,tt] x 
[— 7r,7r]. The gap in this sector is then obtained by 
minimizing Ef roc which yields A = 1 — 2(J X + J y ). 
Note that, in this case, the perturbative result at order 
1 coincides with the nonperturbative result obtained by 
Kitaev' (see also Sec. VIII) and one recovers the transi- 
tion point at J x + J y = 1/2 = J z . 

Results for other sectors (vortex-full, or one vortex ev- 
ery two plaquettes) can also be obtained. We mainly used 
them to check the validity of the coefficients we computed 
perturbatively, as explained in Sec. VIII. 

As a final remark about the 1-QP subspace, let us 
mention a difference with what is obtained when us- 
ing exact fermionization methods. With these methods, 
the low-energy subspace already contains many fermions, 
and one then considers fermionic excitations on top of 
this complicated vacuum to reach high-energy states. In 
our approach, the low-energy states are really empty of 
fermions, and the excitations arc only made of one parti- 
cle, which can thus be qualified of Landau quasiparticlc. 



VII. EFFECTIVE HAMILTONIAN IN THE 
(q > 2)-QP SUBSPACE 

Let us now turn to multi-particle states with the aim 
of showing how the Fermi statistics can be recovered 
from hardcore bosons with a string of spin and phase 
operators. We shall not give many details here, since 
our approach becomes cumbersome when studying multi- 
particle states, and because one knows from exact solu- 
tions that one has to recover free fermions. 



A. Phase factors 

To obtain Fermi statistics, the phase factors appear- 
ing in the string operators S [see Eq. (18)] are of utmost 
importance. These factors do neither appear at the end 



sites of the string operators S [so in particular not at 
all for nearest-neighbor hoppings which are simply the 
t\ operators of Eq. (13)], nor at points of turning back, 
but only at intermediate sites between two truly differ- 
ent other sites. The six possibilities are shown in Fig. 17, 
where the phase factors occur only for the three topmost 
hoppings [phase factors for hoppings not represented in 
the figure can be inferred from hermitian conjugation and 
Eq. (28)]. In the figure, oriented thick (cyan) lines rep- 
resent the string operators S", and the sites j 2 marked 
with a large dot are the ones involving a phase factor 

(— l) b i2 b 32 . As an example, the string operator associ- 
ated to a three-site hopping as the one shown top left in 
Fig. 17 reads 



5'. 



3,3+ni j+2n 1 



In fact, in S.. 



,3l;32,33, 



a phase factor appears at the 



intermediate site i if T. and T. commute, and does 

J 1 3l 32 ' 

not appear if they anticommute. 



B. Fermionic creation operators 

Rigorously, it is impossible to introduce cre- 
ation/annihilation operators for single fermions, because 
fermions should always be created/annihilated in pairs. 
In fact, after choosing a site O as an origin, and after 
choosing a reference path from site O to site i (as was 
done in Sec. VI), the operator (running on this reference 
path) c\ ~ Sq 1 jpl can be considered as a fermionic cre- 
ation operator at site i, once the origin O has been sent 
to infinity (using the same trick as when constructing a 
Dirac monopole in electrodynamics). It should be clear 
from arguments similar to those of Sec. VI that such an 
operator creates a high-energy (spinless) fermion at site 
i, but also creates (or destroys if there is already one) one 
low-energy fermion made of two vortices, top and right 
of site O, as in Fig. 14. It however commutes with all 
other W p operators, except with these two. 

The fermionic anticommutation relations between 
fermion operators at sites i and j can be checked by 
exhausting all possible crossings of two reference paths 
O, . . . ,i and O, . . . ,j. 



C. Multi-particle basis and effective Hamiltonian 

From there on, one can construct a multi-particle basis 
of the Fock-space, as was done for the one-particle basis, 
by successively creating fermions at some sites (after hav- 
ing decided for an ordering of these sites). 

It can then be shown, as was done in the 1-QP sub- 
space, that the Hamiltonian is nothing but a hopping 
Hamiltonian of free fermions in a magnetic field, whose 
flux per plaqucttc is zero or half the flux quantum 
(w„ ~ il)- The phase factors, apart from ensuring 



14 




fkf 




***£ 



FIG. 18: (color online). Illustration of the exchange of two 
particles discussed in the text, for j = i + ri2, k = i + n\ and 
I = i — ri2. 



proper Fermi statistics, also yield the correct expressions 
for the Wp's or product of W p 's, which involve both r's 
and phase factors, and which appear for hoppings around 
closed paths. 



sectors. An alternative route 9,11 ' 12 consists in using the 
Jordan- Wigncr transformation which maps the problem 
onto free spinlcss fcrmions with p-wave pairing. However, 
in both approaches and as is often the case, only periodic 
configurations allow one to obtain analytical expressions 
of the spectrum. In the following, we use Kitaev's ap- 
proach (Majorana fermions) to compute the spectrum in 
several simple periodic configurations characterized by a 

filling factor V = Number of vortex 
° Number or plaquette 

Actually, diagonalizing the Majorana fermion Hamilto- 
nian on this honeycomb lattice ' is completely equivalent 
to analyzing the problem of a free particle on this lattice 
in a transverse magnetic field 45 with a flux per plaquette 
which can take only two values corresponding to w p = ±1 
(see Appendix F for details). The ground-state energy is 
then obtained by filling all levels with negative energy 
which amounts, in a bipartite lattice for which the spec- 
trum is symmetric, to consider half-filling. 

For the three cases considered here, we compute the 
exact spectrum (still assuming J z > J x , J y > 0). Then, 
we perform the perturbativc expansion of the ground- 
state energy up to order 10. This provides some simple 
checks of the results given in Sections V. 



D. An alternative picture 

As was suggested by Levin and Wen in Ref. 42, the 
statistics of the effective quasiparticles can be probed 
with a simple argument. It relies on exchanging two of 
these quasiparticles by using hoppings from the Hamilto- 
nian only, and doing so in such a way that a hopping on 
a bond between two sites as occured exactly once in each 
direction, in order to capture phases coming from the 
statistics only (and not, e. g., from a magnetic Aharonov- 
Bohm-likc phase). 

Let us thus consider the exchange process of two par- 
ticles initially sitting at sites j and I (no other particle 
is present), as depicted in Fig. 18 and whose correspond- 
ing operator sequence is (using only hopping operators 
arising at lowest order) 



-1, 



(38) 



or, equivalently, t % jt\t\ = — tjij^. The sign in the lat- 
ter identity confirms that the quasiparticles made of a 
hardcore boson and an effective spin- 1/2 obey fermionic 
statistics. 



VIII. SIMPLE CHECKS FROM SIMPLE 
VORTEX CONFIGURATIONS 

As shown by Kitaev 7 , the spectrum of the Hamilto- 
nian (1) can be computed exactly by mapping the spin 
system onto free Majorana fcrmions. The main draw- 
back of this mapping is that one has first to work in a 
fixed vortex-sector and, in a second step, perform the 
symmctrization procedure involving all equivalent gauge 



A. Vortex-free configuration v = 

This configuration defined by w p = +1 for all p's is 
of special interest since, in the thermodynamical limit, 
the ground state of H lies in this sector. This is a direct 
consequence of Lieb's theorem for flux phases 41 . The 
spectrum, in this sector, is simply obtained since it is 
equivalent to compute the spectrum of a free particle in 
zero field. The system being periodic with 2 sites per 
unit cell (see Fig. 1), the single-particle spectrum con- 
sists of two bands given by the roots of the following 
characteristic polynomial 



P"=°(e) = e 2 - /(q) 2 , 
where for all q in the reciprocal lattice, 



(39) 



/(q) 2 = 4 J 2 + J 2 + J 2 



J x J y cos (q.(ni - n 2 ) 



J y J z cos(q.n 2 ) + J X J Z cos(q.ni) 



(40) 



The ground-state energy per plaquette is thus given, in 
the thermodynamical limit, by 



^L dqx L d% l/(q)| - 



(41) 



As already found by Kitaev, at the isotropic point 
J x = J y = Jz = 1, one has e v Q =° ~ -1.5746. 

The gap is given by the minimum, in modulus, of 
P" =0 's roots, i. e., min q |/(q)|. Thus, one obtains 



A 



v=0 



2(J ^ Jx Jy) 



(42) 
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Setting J z = 1/2, and considering the perturbative 
limit J z 3> J xl J y , one obtains the following expansion 
for the ground-state energy at order 10 



O 74 

e ° ~ 2 ~T 



5 J 6 



875J 8 3087J 1 



G4 



32 



(43) 



For simplicity, we have set here J x = J y = J . This result 
can be easily recovered by setting w p = +1 for all p's in 
Eq. (21) and using the coefficients given in Appendix C. 

One can also check directly the one-particle spectrum 
by expanding /(q) in the same limit and by comparing it 
with the one-particle spectrum in the vortex-free sector 
obtained in Sec. VI. 



B. Vortex- full configuration v = 1 



The vortex-full sector is defined by w p 



-1 for all 



p's. In the "particle in a field" language, this problem 
corresponds to a magnetic flux per plaquette which is 
half a flux quantum. With the gauge choice shown in 
Fig. 19, the system is periodic with 4 sites per unit cell. 
The single-particle spectrum thus consists of four bands 
given by the roots of the characteristic polynomial 

P v=1 (e) = £ 4 - 8e 2 (4 + Jy + J 2 Z ) + 16. 9 (q) 2 , (44) 

where for all q in the reciprocal lattice, 

.g(q) 2 = 4 + 4 + 4 -2{j 2 J y 2 cos(2q.n 1 )+ (45) 

44cos[q.(m -n 2 )] - 44cos[q.(ni + n 2 )]j. 

The vectors ri\ = (1,0) and n 2 = (0,1) are defined in 
Fig. 19. The ground-state energy per plaquette is given, 
in the thcrmodynamical limit, by 



/ dq x f dq y Jj? c +Jl + j2 

J —IT J — 7T 



1 



l<7(q)l, 



(46) 

Once again, for J x = J y = J z = 1, this expression gives 
eg = 1 ~ —1.5077 in agreement with Kitaev's results' . 

The gap is again given by the minimum, in modulus, 
ofP" 



s roots 



A" =i = 2[ J z -,/J2 + 4 



(47) 



in agreement with results given in Ref. 8. 

As for v — 0, setting J z = 1/2, and considering the 
perturbative limit J z 3> J x , J y , one obtains the following 
expansion for the ground-state energy at order 10 



4 



3J 6 



149J 8 547J 1 



G4 



32 



(48) 



For simplicity, we have also set here J x = J y = J . This 
result can be easily recovered by setting w p = — 1 for all 
p in Eq. (21) and using the coefficients given in Appendix 
C. 
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FIG. 19: (color online). A possible gauge choice realizing the 
vortex-full lattice v = 1. The thin (bold) links are associ- 



ated to Ujk = +1 (Ujk 



T) where j belongs to the black 



sublattice and k to the white one (see Appendix F). The 
eigenvalue of the plaquette operator is then simply given by 



C. Vortex-half configuration v = 1/2 



Let us now consider the vortex-half configuration 
shown in Fig. 20 which is made of alternating vortex- 
free and vortex-full rows. With the gauge choice shown 
in this figure, the system is periodic with 8 sites per unit 
cell. The 8 bands of the single-particle spectrum are given 
from the roots of the following characteristic polynomial 



P"=V2( £ ) = £ 8 _ 16e 6 (J 2 + J2 + j2) + 32e 4[ 3 (4 + J$ + J*) + 4(44 + J V J " + J l J D ~ 2J " J V ^(q.U,)] - 

256 e 2 [4 + 4 + 4 + (4 + 4X4 + 4X4 + 4) - 244(4 + 4) cos^m)] + 

256 14 + 4 + 4 + 444 - 4(44 + J * J y) c° s (q- n i) + 2 4 [44 cos(2q.m) - 44 cos(2q.n 2 )] + 
24{4cos[q.( ni + 2n 2 )] + 4cos[q.(n 1 - 2n 2 )]}|. 
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FIG. 20: (color online). A possible gauge choice realizing 
the vortex-half lattice v = 1/2. Notations are the same as 
in Fig. 19. In this configuration, vortices are localized, in 
alternance on horizontal rows. 
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FIG. 21: (color online). Another possible gauge choice real- 
izing the vortex-half lattice v = 1/2. Notations are the same 
as in Fig. 19. In this configuration, vortices are localized, in 
alternance, on diagonal bands. 



The vectors n\ = (1,0) and n 2 = (0,1) are defined in 
Fig. 20. Note that since the hexagonal lattice is bipartite, 
the single-particle spectrum is even and, consequently, all 
characteristic polynomials are functions of e 2 . Thus, even 
in this vortex-half configuration, one can get analytical 
expressions for the 8 bands since, practically, one only 
has to find the roots of a fourth-order polynomial. 

At the isotropic point, one obtains the ground-state 
energy per plaquette eQ -1 ^ 2 ~ —1.5227. The gap is given 
by the minimum, in modulus, of P"=v 2 'g roots 



A v = 1/2 = 2 j 7 -J.n 



J 2 



(49) 



It is worth noting that the gap, in this sector is exactly 
the same as the one in the vortex- full sector A^ =1 [sec 
Eq. (47)]. 

Expanding the negative roots of P"- 1 / 2 at order 10 
and integrating them out as in the previous sector, one 
gets for J x = J y = J 



v=l/2 
'0 



2 J 4 



3J 6 



A11.P 211J 1 



G4 



32 



(50) 



Finally, one may also consider another vortcx-half con- 
figuration rotated as shown in Fig. 21. The corresponding 
characteristic polynomial is straightforwardly obtained 
from p v=1 / 2 by the permutation J x — ► J y , J y — ► J z , 
J z —>■ J x . However, since the perturbation is performed 
in the limit J z 3> J x , Jy it leads to a different expression 
for the expanded ground-state energy. In this case, one 
gets for J x = J y = J 



v=l/2 _ _ 1 _ j2 _ J 4 _ J S 

~ 2 4 4 



109J S 
64 



59J 



10 



16 



(51) 



Once again, both expressions (50) and (51) can be re- 
covered from Eq. (21) using the coefficients given in Ap- 
pendix C. 

The various results obtained for v = 0, 1, 1/2 provide 
(partial) checks of the coefficients given in Appendices C 



and D and show the power of the PCUTs to compute 
high-order expansion for the spectrum. In the next sec- 
tion, we shall show that this method is also an efficient 
tool to tackle more complex problematics. 



IX. OBSERVABLES 

One of the advantages of the CUTs method is that 
it allows one to obtain the effective form of any obscrv- 
ables and to compute its matrix elements in the eigen- 
basis of the Hamiltonian. The aim of this section is two- 
fold. First, we compute pcrturbatively the spin-spin cor- 
relations and show that they admit a plaquette-operator 
expansion similar to that of the spectrum. The second 
part of this section is dedicated to the most fundamental 
problem of local spin operations onto the ground state. 
Following Ref. 16, we show that single-spin operations 
create anyons but also fcrmions. We compute the spec- 
tral weights of various states stemming from such opera- 
tions and we also analyze the action of string operations 
which allow for manipulation of anyons. Finally, we give 
a procedure to derive the operators which create anyons 
without fcrmions and show that they involve tricky su- 
perpositions of multi-spin operators. 



A. Spin-spin correlation functions 

The Hamiltonian (1) is invariant under the time- 
reversal symmetry since it is a quadratic function of the 
spin operators. Thus, any expectation value of an odd 
number of spin operators vanishes (such as the mag- 
netization (erf)). Note also that the absence of odd 
cycles ensures that the eigenstates do not break this 
symmetry 7,26 ' 27 . 

In addition, the only nonvanishing spin correlators are 
those involving products of af on a-dimers 12 ' 17 ' 18 . In 
this section, we focus on the spin-spin correlation func- 
tions and their expression in the 0-QP sector. More pre- 
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cisely, we consider the following operators = crfa^ 
where (i,j) is an a-dimer. To compute these quantities, 
we proceed in a way similar to what we have already done 
to derive the effective Hamiltonian : 

• we express the observable in the ESB language ; 

• we compute its effective form perturbatively as ex- 
plained in Appendix B (see also Ref. 45 for a de- 
tailed discussion) ; 

• we project it out in the sector of interest. 

Using the ESB form of the spin operators (8), we 
straightforwardly achieve the first step mentioned above 
for the three correlation functions 

^i>iW,. = (&i + &iK? +ni (&t +m +&. +m ), (52) 

<o<. = l-2&t 6l = (_i)*K. (54) 

To avoid any ambiguity, we keep track of the type of 
sites (• or o) but we are working, at this stage, on the 
effective square lattice. Next, we turn to the second step 
using the perturbative expansion described in Appendix 
B. In the present case, we pushed the calculation up to 
order 6 and, finally, we focus on the 0-QP sector. 

As one expects, the effective form of the spin-spin cor- 
relation function is similar to that of the effective Hamil- 
tonian. This is due to the fact that, in the low-energy 
sector, W p s are the only degrees of freedom. Thus, we 
obtain, an expansion in terms of the plaqucttc operators 

C%*\ q = = a aa - ]T b^_ p W Pl W P2 ...W Pn . 
{pi,...,p„} 

(55) 

The coefficients a aa and b p " Pn are given in Appendix 
E up to order 6. Here again, we can see that these corre- 
lation functions involve interactions between connected 
or disconnected plaquettes. 

As a simple check of our expression, one can eas- 
ily compute the nonperturbative correlation function in 
the vortex-free and the vortex-full sector thanks to the 
Hellman-Feynman theorem. Indeed, in these sectors, all 
sites are equivalent so that one readily gets the expecta- 
tion value 

r)p v 

(C?f*\ g=0 )„ = C?f<\ v q=0 = - I £. (56) 

for both cases v — 0, 1 for which the ground-state ener- 
gies are given in Eqs. (41) and (46). As in the previous 
section, the subscript v indicates that we consider the 
ground state of the sector with filling factor v = 0,1. 
Then, expanding these expressions (before derivation), 
setting J z = 1/2 and for simplicity J x = J y = </, one 
gets 

q 74 

CIXZI = 1-2J 2 - — -25 J 6 , (57) 
Ct*\ v =l = J + -f- + if - = C^=l (58) 
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FIG. 22: (color online). Behavior of the spectral weights 1^ 
for fermion numbers n = 0, 2, 4, as a function of the coupling 
J — J x = J y for J z — 1/2. Gray plaquettes in the insets show 
the positions pf and p| at which the anyons are created under 
the action of r* . 

for the vortex-free sector and 

CtXzl = 1 - 2 J 2 + — - 15 J 6 , (59) 

CtKZl = J Y + -f- - Cf! |^ \ (60) 

for the vortex-full sector. As can be checked, these results 
can be recovered using the coefficients given in Appendix 
E and Eq. (55). We emphasize that, as for the spectrum, 
our expressions allow us to investigate arbitrary vortex 
configurations such as sparse vortex ones, recently stud- 
ied numerically 1 ' 23 . 



B. Creation of anyons 

Let us now analyze the action of a single-spin opera- 
tion onto the ground state and following Ref. 16, let us 
focus on a? As for the correlation functions, one first 
has to write this operator in the ESB formalism which 
is, again, straighforward since erf. = r?. Then, one com- 
putes its renormalization under the unitary transforma- 
tion U which "diagonalizes" the Hamiltonian. Finally, 
one can compute any matrix element of this observable 
between any eigenstatcs. 

At order 0, the observable is not renormalized and one 
has U>t?U = t\. When this operator acts onto the 
ground state which is in the vortex-free sector, it thus 
simply flips the two plaquettes as shown in Fig. 22. In 
other words, it creates two anyons and nothing else. 

At order 1, one gets 

C/V? U = r? [1 + (J x vt_ ni + J v vt_ n2 + h.c.)] , (61) 
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showing that things arc more complex since pairs of parti- 
cles (fermions) are created. It means that, at this order, 
t? couples the 0-QP subspace of the vortex-free sector 
with the 2-QP subspace of the two-vortex sector dis- 
cussed above. To have a physical quantitative picture 
of such processes, let us compute the spectral weights 
defined as 



V 



k 



\{{pl,pl),n,k\ 



(62) 



where \{p}, n, k) denotes the eigenstate of H in a sector 
given by an anyon configuration w p = — 1, and n high- 
energy quasiparticles with quantum numbers k. Here, 
the plaquettes p\ and p\ are as indicated in the inset 
of Fig. 22. This quantity measures the weight of all n- 
fermion contributions obtained by the action of r? onto 
the ground state |0) which contains no fermion and no 
anyon. As it should, these spectral weights satisfy the 
1. At order 6, one gets 



sum rule V n 
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(63) 



(64) 



which shows the importance of the two-fermion states 
for increasing couplings as can be seen in Fig. 22. Note 
that the sum rule is fulfilled here implying i^ >4 = at 
order 6. Actually, one may consider representative curves 
in Fig. 22 as almost converged since order 8 corrections 
would bring very small corrections. 

To summarize, one must realize that local spin oper- 
ations onto the ground state create anyons (here two) 
but also give rise to fermionic excitations whose weight 
increases significantly with the perturbation. 



C. Manipulation of anyons 

Another important question concerns the manipula- 
tion of the anyons which, as shown above, may be created 
by local spin operations. Such an issue is of special in- 
terest for experiments aiming at braiding anyons . This 
topic has been the subject of a recent controverse with 
Zhang et al. 33 ' 46 who completely neglected the existence 
of fermions in this model. Following Jiang et al. 30 who 
proposed an ingenious protocol to detect anyons statis- 
tics, we wish to compute the action of a string operator 
onto the ground state. 

For simplicity, we consider here the operator 
S = ria=i m a i • a l° n g a horizontal line of the original 
brick- wall lattice (see Fig. 23 with to = 3 for notations). 

At order 0, it is simple to see that S first creates two 
anyons and make one of them jump in the direction of the 



FIG. 23: (color online). Action, in the vortex-free sector, of 
the string operator S = cr^.of^.of^.. Each operator flips 
the two plaquettes adjacent to the z-dimer it is attached to 
but also creates fermionic excitations (not shown). 



string so that, at the end, one eventually has one anyon 
at plaquette 1, another anyon at the plaquette m + 1, 
and no fermion. 

However, at higher orders, as previously, such an op- 
eration creates fermions. To quantify this phenomenon, 
we consider the probability V = |({l,m+ 1},0|S'|0)| 2 to 
find the final state in the lowest-energy state (no fermion) 
with anyons at plaquettes 1 and (to + 1) which coincides 
with I§ for jn = l. In Rcf. 16, we computed this prob- 
ability at order 2, but here, we go beyond and give the 
result at order 6 



V 



1 - m( j 2 + J 2 ) + (65) 

{J$ + J 4 ) + (to 2 - 8to + 3) J 2 X J 2 - 
to(to 2 - 12to + 32) 
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The main reason to perform this high-order calcula- 
tion is that the above expression pleads in favor of an 
exponentiated form linear with to. Indeed, although we 
have no proof, we conjecture that V can be recast into 
exp(A — to-B) as suggested in footnote 4 in Ref.' . It is 
indeed striking to see that the expression (65) which is a 
polynomial of the variable to, can be seen as the expan- 
sion of such a simple form with, at order 6 



.4 
B 



iJ X Jy 

J x + jy 

n44 



2Q{JU : 

+ 8J 2 j'; 



J-x Jy ) 



- J x j y ), 

2(4 + Jy) + 

■^(4 + 4)- 



(66) 



(67) 



Further, it is clear that V is bounded by and 1 for any 
to, which is clearly not the case if one considers (65). 
Let us also note that the fact that U't?U is found to 
be proportional to r? [see Eq. (61)] strengthen the idea 
of an exponential form of this effective observable and 
hence for S. 

We display the results at various order in Fig. 24 us- 
ing the expanded form (65) and the exponential form. 
As can be clearly seen, the exponential form seems to 
be well-behaved. In addition, the order 6 expansion of 
A and B seems to provide an almost converged result 
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FIG. 24: (color online) 
computed for m — 25. Top : nonresummed (bare) expression 
(65) at order 2 [red (bottom)], 4 [green (top)] and 6 [blue 
(middle)]. Bottom : exponentiated form exp(A — mB) at 
order 2 [red (top)], 4 [green (middle)] and 6 [blue (bottom)]. 



when put in the exponential. Thus, we claim that one 
can use this form to obtain a very accur 



experiments 
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D. Anyons without fermions 

As discussed previously, local or string spin operations 
create fermions. However, in experiments, one wishes 
to manipulate anyons without being spoiled by these 
fermions . In other words, the ideal operations would 
consist in exciting plaquettes (only) while remaining in 
the ground state of the corresponding vortex sector. In 
this section, we shall show, perturbatively, that it is pos- 
sible to do so even if the form of such operators is hard 
to implement in realistic devices. 

As an example, let us determine the operator f2^ cre- 
ating two vortices at the left and right plaquettes of a 
given site i [see inset Fig. 22 (left)]. This operator must 
be such that fli c g = U'CliU = r? which indeed leads to 
Iq = 1 . Note that this procedure is the inverse of what is 
usually done with CUTs since, here, we wish to compute 
the bare observable given the effective observable instead 
of the opposite. 

Let us assume that this operator has a perturbative 



expansion, namely 



where contains all operators of order k and thus 
associated to J l x J™ (with I + m = k). At order 0, op- 
erators are not renormalized so that one obviously has 
= t? . The rcnormalization of fij under the unitary 
transformation U reads 

n icS = £ nfl = £ L/tnf ) u = £ £ n$-M , (69) 



(68) 



where Q 



he: 



(k),[l) 

l eff 



A' ei 



is of order (k + I). Since, at order 0, one 
has fiioff = t?, one must have, at each order r > 

EE^=0 (70) 



where the sum is restricted to values of indices such that 
k + I = r. At order 1, this leads to 



(0),[i] 

l eff 



n,^ = o. 



(71) 



Using Eq. (61) and the fact that f^' 1 " 1 = one 
then obtains 

«f 5 = - (J*vt- ni + Jy4- n2 + h.c) . (72) 
Using the inverse mapping of Eq. (8) 

(73) 
(74) 
(75) 

(76) 

one finally gets, in the original spin language and at order 
1, 



value of 
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= a i.» (j i, 


braiding 
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1 (rr* 
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Jx((T. 



1 r 
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(77) 
(78) 



T I Z X X 

J V \ a i-n2,» a i-ri2,o a i,' 



V V z 

°i-n 2 ,o0i,.^i,o 



This expression shows that to create anyons without 
fermions, one has to build a complex superposition of 
operators with fine-tuned coefficients. At order 1 consid- 
ered here, such states require single and triple spin-flip 
operations but, of course, higher-order corrections would 
involve higher order spin-flip processes. Such constraints 
makes creation of anyons without fermions via local op- 
erations difficult experimentally . 
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X. CONCLUSION AND PERSPECTIVES 

We have analyzed perturbatively the gapped phase of 
the Kitaev honeycomb model in the isolatcd-dimer limit 
using the continuous unitary tranformations. We have 
thus derived the low-energy effective theory up to order 
10 which has been found to describe an interacting anyon 
system. This result has to be contrasted with the order 
4 result which predicts a free anyon system'. We also 
showed that the excitations in each vortex sector obey 
fermionic statistics. 

In a second step, we focused on the action of local 
spin operators onto the ground state and we have shown 
that they generate both anyons and fermions. We also 
gave the form of the operator which creates anyons with- 
out fermions. This operator involves multi-spin operators 
which may be hard to implement experimentally. 

Of course, several questions remain open in this model. 
As explained by Kitaev, there exists a gaplcss phase 
which is associated to non-Abelian anyons. The influ- 
ence of a magnetic field in this phase is certainly one of 
the most challenging question and should reveal rich phe- 
nomena. Note that the effect of a magnetic field in the 
toric code already gives rise to a nontrivial phase diagram 
as recently discussed in Ref. 47,48. 

Another interesting issue concerns the time evolution 
of local excitations. Indeed, in Sec. IXC, we discussed 
the effect of a string operator onto the ground state but 
we always considered static quantities. Although exper- 
imentally, succesive spin operations may be performed 
on "short" time scales, it would be of primer interest 
to compute the spreading of fermionic excitations during 
the braiding processes proposed to detect anyons ' . 
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APPENDIX A: STRUCTURE OF THE 
EFFECTIVE HAMILTONIAN 

As we have seen in Sec. Ill, when setting J z = 1/2, 
the Hamiltonian (1) can be written as 

N 

H = --+Q + n + T +2 + T_ 2} (Al) 

where N is the number of z-dimers, Q is the particle- 
number operator, To contains the pure hopping opera- 
tors which does not change the number of particles and, 
T + 2 (T-2) creates (annihilates) pairs of particles. The 
operators Tq ±2 are proportional to the small parameters 
from which the perturbation theory is performed. 

The idea of the present approach is to transform the 
Hamiltonian (9) into an effective one which conserves 



the particle number. Of course, in general, this cannot 
be achieved exactly and, as often, one has to perform 
a perturbative expansion. To achieve this goal, a very 
powerful tool is the continuous unitary transformations 
method 34 . For the problem at hand, Knetter and Uhrig ::l 
have developed a code which computes the coefficients of 
this expansion at high orders 11 . Practically, one must 
keep in mind that, at order 10 which is the maximum or- 
der considered in this paper, one already has more than 
10 4 terms. We refer the interested reader to Ref. 39 for 
a detailed derivation and we give below, for illustration, 
the results up to order 4. 



Order 


Operator O 


Coefficient c 




1 


To 


1 


1 


2 


T-2T+2 


-1/2 





2 


T+2T-2 


1/2 
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3 


T-2T0T+2 


1/4 





3 


ToT— 2T+2 


-1/8 


1 


3 


T-2T+2T0 


-1/8 


1 


3 


T+2T-2T0 


-1/8 


2 


3 


T T +2 T-2 


-1/8 


2 


3 


T+2T0T—2 


1/4 


3 


4 


1 0/ o/io/io 
1 — A — Z -*- + Z -*- -\- Z 


-1/16 





4 


T-2T0T0T+2 


-1/8 





4 


T-2T+2T-2T+2 


1/8 





4 


To T_ 2 To T+2 


1/8 


1 


4 


To To T- 2 T+2 


-1/32 


1 


4 


T-2T0T+2T0 


1/8 


1 


4 


To T_ 2 T+2 To 


-1/16 


1 


4 


T_ 2 T+2 To To 


-1/32 


1 


4 


T+2T-2T0T0 


1/32 


2 


4 


To T+2 T_ 2 To 


1/16 


2 


4 


ToT)T+2T_2 


1/32 


2 


4 


T+2T_2T+2T_2 


-1/8 


2 


4 


T+2T0T-2T0 


-1/8 


3 


4 


To T+2 To T_ 2 


-1/8 


3 


4 


T+2T0T0T-2 


1/8 


3 


4 


T+2T+2T-2T-2 


1/16 


4 



TABLE I: Operators appearing in H e a with its corresponding 
coefficient up to order 4, together with the g-particle subspace 
they start to act on. 

The operators and corresponding coefficients are put in 
Table I, together with the lowest number g m j n of parti- 
cles such that the operator has, a priori, a nonzero action 
within the g-particle subspace for q q m in- 9min is found 
by requiring that the number of particles in the system 
is always positive, and by using the fact that T projects 
out 0-particle states. Note that some terms may vanish 
for more subtle reasons. For example, the third order 
term T_2ToT+2 does not act on the 0-QP states. Indeed, 
T+2 creates a 2-QP state ; then Tq makes one of the 



particle hop ; and finally T_2 tries to annihilate two par- 
ticles, but cannot, since these are not nearest-neighbor 
anymore, due to the hopping. 

One can then directly write the effective Hamiltonian 

H eS = -— + Q + ^2c i O i , (A2) 

i 

where Oi is the i th element of the column "operator" of 
Table 1 and c; the associated coefficient, the order being 
given by the first column. By construction, the effective 
Hamiltonian conserves the particle number and the en- 
ergy states are ordered according to their quasiparticle 
number, the ground state being in the 0-QP sector. Fur- 
thermore, since [H s, Q] = 0, one may also rewrite the 
effective Hamiltonian in the following form 

H e s = ^-ffofflg, (A3) 

where i? e ff|g denotes the projection of H c h onto the 
q— QP sector. Note that it is not the decoupling used 
in the CUTs community where usually one gathers all 
operators which contain exactly q creation and q annihi- 
lation operators and thus act on the q'-QP sector with 

q'>q- 

Finally, one must analyze each sector defined by the 
number of quasiparticles and determine the action of each 
operator Oi in the corresponding subspace. This is the 
nontrivial part of the job which depends on the problem 
under consideration. Let us emphasize that if each oper- 
ator only starts to act in the <7min-QP sector, it has also, 
in general, a nontrivial action on the q-QP sectors for 
q > Qmin- 

APPENDIX B: PERTURB ATIVE EXPANSION 
OF OBSERVABLES 

In this appendix, we give the general perturbative ex- 
pansion of any observable fl obtained with the CUTs 
using the quasiparticle number conserving generator. As 
is the case for the Kitaev model, we suppose that the 
Hamiltonian of the system can be casted in the following 
form 

H = Q + T- 2 +T Q + T+ 2 , (Bl) 

and satifies the hypothesis given after Eq. (16). In this 
case, the flow equations obtained from the CUTs method 
can be solved perturbatively 4 ', and the effective observ- 
able can be written as: 

n eS = n + Y^aOi, (B2) 

i 

where Oi is the i th element of the column "operator" of 
Table II and c, the associated coefficient, the order being 
given by the first column. 
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Order 


Operator 


Coefficient 
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T- 2 fi 
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t+2 n t+2 


-1/4 


2 
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1/8 
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VT- 2 T 
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2 
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-1/8 
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T T-2 


1/4 


2 
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1/4 


2 


^ T+2 T_2 


-1/8 


2 


n T+2 To 


-1/4 


2 


n t+2 t+2 


1/8 



TABLE II: Operators appearing in Q e ff with the correspond- 
ing coefficient up to order 2. 

At order 6 considered in this paper for the correlation 
functions, there are several thousands of terms to con- 
sider. Once this effective form is derived, one then has to 
analyze it in the quasiparticle sector of interest as done 
for the effective Hamiltonian. 



APPENDIX C: COEFFICIENTS OF THE 
PERTURBATIVE EXPANSION OF THE 
HAMILTONIAN IN THE O-QP SECTOR 

As explained in Sec. V, the effective Hamiltonian in 
the 0-QP sector schematically reads 

H e s\q=Q = Eq - Cpi,...,p„Wpi Wp2 . . . W pn . 

{pi,...,p»} 

(CI) 

where {pi,p%, . . . ,p n } denotes a set of n plaqucttes and 
W p are conserved plaquettc operators. The form of 
the effective Hamiltonian is translationaly invariant (of 
course the configuration of the Wp's need not be !), so 
that C pi) . Pn in fact only depends on relative coordi- 
nates of the plaquettes, and we will use (except for the 
one-plaquette coefficient) the notation Cp 2 - Pl ^.^ Pn - Pl = 
Cpi,...,p„- Here, we give the perturbative expansion up 



to order 10 of Eq and the C's in the limiting case Three-plaquette terms 

Jx,J y "C J z - Setting J z = 1/2, one gets the following 

results. 
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where TV is the number of z-dimers. 
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Two-plaquette terms 
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Four-plaquette terms 
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As can be seen from these expansions, the number of 
interacting plaqucttes increases with the order of the per- 
turbation theory. 



APPENDIX D: COEFFICIENTS OF THE 
PERTURB ATIVE EXPANSION OF THE 
HAMILTONIAN IN THE 1-QP SECTOR 
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Five-plaquette terms 



In the 1-QP sector, the effective Hamiltonian reads [see 
Eqs. (26)-(30)] 



H c s\ q= i = H c s \g—o + n — ^ Dj_ Lt 



■ t 3n 



3i 



(Dl) 

where the sum is performed over all non self-retracing 
paths of length n starting at site j 1 and ending at site 
j n . The operators t\ arc defined in Eqs. (10-13). Since 
the D's do not depend on the initial site, we introduce 
Dj —j , t j —j = Dj it ___j . From the symmetries of 
the underlying lattice, it is clear that we can limit the 
analysis to processes involving a first jump in the direc- 
tion +m or +TI2. 

We give below the perturbative expansion of /x and 
the -D's in the limiting case J Xl J y <C J z up to order 4 
and set J z = 1/2. Note that one could reach order 10 
as for the 0-QP sector if needed. However, as explained 
in Sec. VI, it is simpler, in this sector, to use directly 
the Majorana formalism which is nonperturbative and 
requires a comparable numerical effort. 
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Chemical potential 



H = 1 + Jl + 



One-hopping terms 



1 



7 7 J 7 j z 



1 

2"» 



n — T — - T 3 J 2 T 



(D2) 



24 



Two-hopping terms 



Four-hopping terms 
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FIG. 25: (color online). Labeling of the plaquettes by a site 
index (i) and a position it, d, Z, r with respect to that site. 



clockwise or anti-clockwise but the product of leads 
exactly to the same operator 6tfe i W p . 
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APPENDIX E: COEFFICIENTS OF THE 
SPIN-SPIN CORRELATION FUNCTION IN THE 
0-QP SECTOR 

As discussed in Sec. IX A, the spin-spin correlation 
functions C"^ = crfa^ computed on any eigenstate of 
77 is nonvanishing only if a = j3 and if i and j belong to 
the same dimer which is of a type. 

We give below the perturbative expansion of these cor- 
relation functions in the 0-QP sector. 
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In the perturbative approach we use, note that a z- 
dimer in the honeycomb lattice becomes a single site i in 
the effective square lattice. 

As for the Hamiltonian in the 0-QP sector [see 
Eq. (21)] , we obtain an expansion which can be expressed 
only in terms of the plaquette operators, namely 



c: 



1 9=0 



V b zz 

/ J Pi, 

{Pi, ■■■,Pn] 



., Pn W pi W P2 ...W Pn . (El) 



Below, we give the results up to order 6 and we index a 
plaquette p by a site i and an indice u,d,l,r according 
to notations given in Fig. 25. 



Constant term 



Additionally, there are some terms corresponding to 
processes where the particle hops one times around a pla- 
quette. Note that the plaquette involved can be covered 
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Coefficients of Cfj 



Contrary to z-dimers, x-dimcrs remain dimcrs pcrtur- 
bativcly. Here again, one obtains an expansion in terms 
of plaquettes for these observables, in the 0-QP sector, 
which can be written as 



°i,i+ni 1 9=0 



E 



b x P t..., p W pl W P2 ...W Pn , 



(E2) 



(E3) 

We give below the expansion of the coefficients up to 
order 5 (only odd orders are nonvanishing) and, as pre- 
viously, we index a plaquette p by a site i and an index 
u, d, I, r according to the notations given in Fig. 25. In 
the following we consider a dimer located at (i, i + rii). 



Constant term 



Two-plaquette terms 
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Two-plaquette terms 
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(E4) 

3. Coefficients of 

The correlation functions are straightforwardly 
obtained from Cfj by exchanging directions a; and y as 
well as J x and </„. 



APPENDIX F: CORRESPONDENCE BETWEEN 
THE MAJORANA FERMION SPECTRUM AND 
A FREE-PARTICLE PROBLEM IN A MAGNETIC 
FIELD 

As shown by Kitaev' , the spin Hamiltonian (1) can be 
mapped onto the following Majorana fermion Hamilto- 
nian 

where A is a skew-symmetric matrix of size 2N x 2N (N 
being the number of plaquctte) and where the Cj's are 
the (hermitian) Majorana operators which obey c| = 1 
and CjCfe ~ —c k Cj if j =/= k. The sum is performed over 
all sites j and k of the honeycomb (brickwall) lattice and 

A jk = 2J a u jk , (F2) 

if the link (j,k) is of a- type and otherwise. The Uj k 's 
are antisymmetric (uj k = —Ukj) and take the values ±1. 
These numbers define the vortex configuration through : 
w p = Y\(j i.) 6 p Ujk where j belongs to the black sublattice 
and k to the white one (see Fig. la). We refer the inter- 
ested reader to Ref. 7 for details. In the very end, the 
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whole spectrum of H can be obtained once one knows 
the spectrum of iA, e. g. the ground-state energy per 
plaquette is given by 

e = -^Tr|U|. (F3) 

For a bipartite lattice such as the honeycomb lattice, 
we shall now show that the spectrum of \A is the same as 
the one-particle spectrum of the following Hamiltonian 

H' = -^J2 A 'jk(a]a k + h.c), (F4) 

where (a •) are standard spinless fermion creation (an- 
nihilation) operators. The Hamiltonian H' describes free 
spinless fcrmions hopping in a honeycomb lattice in a 
magnetic field with a flux per plaquctte which equals 
zero (w p = +1) or half a flux quantum (w p = — 1). The 
spectra of fT (with one fermion) and of iA are identical 
provided 

A' jk = 2J a u' jk , (F5) 

with u'j k = +u' k j- The choice of the u'j k is as previously 
dictated by the flux configuration via w p = Ylu k ) ep u 'jk 
where, in this case, the u'j k are not oriented but still take 
the value ±1. 

To show this, consider an eigenstate \ip) of the ma- 
trix iA with energy E and let us denote tpj = (jIV') hs 
component on site j. This state satisfies 

^2\A dk il> k = Ei/> j , (F6) 
k 

Since the honeycomb lattice is bipartite, we can always 
set Uj k — u'j k if j is a white site (and k black), and 
u jk = ~ u 'jk ^ j i s a black site (and k white). Then, one 
can easily check that the state \<j>) defined by <f>j = —ipj 
if j is a black site and <j>j = — i ipj if it is a white site, 
satisfies 

-J2 A ' j k<l>k=Ecf> j , (F7) 
k 

so that \(j>) is an eigenstate of H' with the energy E. This 
shows that H' (with one particle) and iA are isospectral. 
We insist on the fact that this correspondence only holds 
for a bipartite lattice but is no longer true in the presence 
of odd cycles. 
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